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ABSTRACT 


This thesis presents a Fortran program that numerically solves the 
steady-state matrix Riccati equation of the quadratic cost optimal 
control problem. Each step of the program is presented, analytically 
and computationally. The check points incorporated in the program and 
the input parameters that can be used to assure a correct solution are 
identified and discussed. Difficulties encountered when verifying the 


program, and the suggested solutions, are also presented. 
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I. INTRODUCTION 


The Fortran computer program presented in this thesis provides a 
relatively quick and reliable means for determining the unique, 
symmetric, steady-state solution P of the nonlinear matrix Riccati 
equation: 


] 


P=0=- PA-A'P - C'C + PBR B'P (1) 


occurring in optimal control theory. The remaining equation variables 
are defined in the following statement of the quadratic cost optimal 
control problem. 

Given the linear, time-invariant, completely controllable system 


defined by the state equations: 


x (t) = Ax(t) + Bu(t) (2) 

x (0) = Xo (3) 
where: 

x(t) is ann x 1 state vector 

u(t) is an mx 1 unconstrained control vector 

A is an n X n matrix | 

B is ann X m matrix, 
determine the control vector u*(-) which minimizes the quadratic cost 


functional: 


CHET Ge Pie f PeGuGlec et ust t)Ru(t) jdt (4) 
QO 


where: 
Cis anm xn matrix 


R is an m xm symmetric, positive definite matrix 


and the matrix pair (A,C) is completely observable. 
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It has been shown [Ref. 1] that u*(-) is given by the linear feed- 
back law: 


1 


itt e—eeh  BUPx(t) = =L*x(t) (5) 


where P is the unique solution to matrix equation (1). 
For the system of equation (2) to be completely controllable, the 


n X mn augmented matrix G, 


G = (B | AB] AcB [+> | ADB), (6) 


must contain n linearly independent column vectors, or equivalently, 
the rank of G must be n. For the matrix pair (A,C) to be completely 


observable, the n x mn augmented matrix H, 


i e@omecrmir )ccimiess | (ar jt let). 
must contain n linearly independent column vectors or the rank of H 
must be n. 

When u*(-) is given by equation (5), equation (2) can be rewritten 


a“ x(t) = Ax(t) - BL*x(t) = (A - BL*)x(t) (7) 


which, using equation (3), has the general solution: 


y= fA - BLA)E (8) 


x(t 
Because the system is completely controllable, as time approaches 
infinity x(t) must remain bounded. To satisfy this requirement, the 
matrix (A - BL*) must be a stable matrix or, alternatively, all the 
eigenvalves of the matrix (A - BL*) must have negative real parts. As 


a consequence of this: 


lim 7 lim (A - BL*) _ 
pees Xa € =a (3) 





To evaluate the quadratic cost associated with the optimal con- 


trol u*(t), substitute equation (5) into equation (4) and write: 


J* = J(x,3u*(t)) = f be cercrente) + x'(t)L*'RL*x(t)] dt. (10) 
0 


Substituting equation (8) for x(t), expanding L*'RL*, and factoring out 


the common terms: 


Je = x! | f{ ofA - BL*)'troic 4 pprlpipy (A - BL*)t gy xo» (11) 
6) 


1 


From equation (1) C'C = -PA -A'P + PBR B'P or: 


1 


CuC sepeR sbeP = -P(A = BL*) - {A = BL*)'P . (12) 


Substituting equation (12) into equation (11) and integrating by parts 


the first term of the resulting integral expression, we get: 


ye “ees r-P(A ~ BL*)] (A - BL*)t 
6) 
- . Q(A - BL*)'t PI a(A - BL*)t - 


; fh (eau) e(A - BL*)'t pp) (A - BL*)t Gy (13) 
O 


Using equation (9) to evaluate the upper limit of the first term, the 


first term reduces to P. Recognizing: 


Zt a Zt 


eee) ght. =e) 2 (14) 
i=0 °° 
we see that the second term in equation (13) will cancel the second 


term in equation (11) when equation (12) is substituted. Thus: 


d= ssa (15) 





IT. ANALYTICAL METHODS OF SOLUTION 


A. ANALYTICAL METHOD OF KLEINMAN 
The main computer program of this thesis is based upon a method for 
solving the steady-state Riccati equation published by David L. Kleinman 
[Ref. 2] in 1968. The iteration scheme for solving euqation (1) is 
presented in the following theorem by Kleinman. | 
1. Kleinman's Theorem 
Pea ea eenmlianc..... be the n x n (unique) positive definite 


matrix solution of the linear algebraic matrix equation: 


0 = AVE + VpAy +C'C + LyRLy (16) 
where, recursively, 

L, = RUBY, k = 1,243)... 

ty = tal = iE KS BOC i os. 


and where L. 1s chosen such that the matrix ah = A - BLY has eigenvalues 


with negative real parts. 


Then: 1) P< View] < Vy ee ere 0, 2 os. 
lim _ 
2) ae ae 


The notation X > Y [X 2 Y] means that the matrix X - Y is 
positive [semi] definite. 
2. The Cost Matrix 
The proof of this theorem requires the introduction and defini- 
tion of a cost matrix Vi: Assume that uy (x(t) ) = - Lx({t) is an arbitrary 


feedback law, with feedback gains of L, and uy (x(t)) is applied to the 





system of equation (2). Following a development similar to equations 


(7), (8) and (15), the resulting quadratic cost function is: 
J(xX3u, (x(t) )) = XV Xo 


where Ve is the cost matrix associated with the feedback gains L and is 


given by: 


v= [ (A - BL)'t rere a Lepiy ofA ~ BLIt ge | (17) 
O 


Vv is bounded if and only if the closed-loop system control matrix (A - BL) 
is stable. If oT is bounded, it becomes the unique (positive definite) 


solution of the linear matrix equation: 


O = (A - BL)'V, + Vi CA - BL) +C'C +L'‘RL . (18) 


L 
Examining the first term of equation (18) and substituting 


equation (17) for Vi» we can verify this relationship: 
(A - BL)'V, = f (A = git et ~ BL) tree ¢ Liptay et - SLIt ay. 
0 


Integrating by parts and using equation (14) we have: 


c 


A - BL)t 


A - BL)'t 1 i ( 
[C'C + L'RLJ e t=0 


(A - BL)'V, = e | 


1 


2 e(A - BL)'troic 4 Lipy @fA - BLIt (a _ Bly at. 
O 


Using equation (9) to evaluate the upper limit of the first term, the 


first term reduces to C'C + L'RL and we have, 
(A - 3109 a eee Oo wie RL = Vi OA - BL) , 


the desired result. 





Recalling the result of equation (15) we see that for the 
optimal control of equation (5) we have: 


ieekh ¢ (19) 


teh Ly and L. are the gain matrices associated with the cost 


matrices V, and V4, it can be shown [Ref. 3] that: 


17 \o af i [(L) - bo) 'R(Ly - Lo) 


- (Ly - Lo)’ (B'V, - RL») 


(ies BL, )t 


(B'V, - RLy)'(L, - L,)] e dt (20) 


or: 


V)- V5 = J ef - BL,)'t [(L; - Lp)'R(Ly = Ly) 


- (Ly) - Ly)'(B'Y, - RLo) 


et.) etsy) e'*  Py)t at. (21) 


If either matrix (A - BL.) or matrix (A - BL) is unstable, Vy 
or Vo» respectively, will be unbounded and care must be exercised in 
using the above formulas. 

3. Proof of Kleinman's Theorem 

1D ee ee Vig] < VE coe kes Os1g2.°°*, Let Ne be the cost 

matrix for Lo: the initial guess that yields a stable closed-loop 


] 


system control matrix, and let Vy be the cost matrix for Ly = Ro BV): 


Substituting V, and V, into equation (20) and noting: 
et = (ek y 
R = Y'Y, where Y is unique 


Baa - RL = 0 





= A,t A,t 


Define Z(t) = Y(L, - Len. It has been shown [Ref. 4] that for 


Z(t) a real matrix: 


(ae es) O fopeotletes (), 


Thus v = Wes f Zt) 7 eiadt <> 0 sor: 
O 


ee que (22) 


Let V_ be the cost matrix associated with L*, use equation (19) 


and substitute into equation (21). Noting that: 


BP = RL* = 0 


V,- P= f aye [(Ly - L*)'R(L, - L*)] pee aye 


lias time define Z(t) = Y(L, - Lt)ehyt and we have: 


] 


ie 2) 


a Pp = f{ ee vatinct > 010K: 


0 
Ue Be (23) 
Combining the results of equations (22) and (23) we get: 


ls 2 Vy paul 


meaning that Vy 1s bounded above by P and below by WAC Thus, A, is a 
Stable matrix and V, Satisfies equation (16) with k = 1. Similar 
arguments can be made for k = 2,3,4,... yielding: 


ve > Vy > Vo eee Vp > Veg =P 


the desired result. 
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2) Vim Vy era ine he VE = V exists by a theorem on monotonic 


convergence of positive operators [Ref. 5]. Thus, in the limit V, = 
] -| 


Vee] and AL =A - BL = A - BR B Vy = A - BR BV E47: Since Ay = 
A' - VBR” |B! equation (16) becomes, in the limit as k approaches 
infinity, 
= ' -lo lee j lag 
O=A VE - VBR B Vp + VLA - VBR B Vv. + CG ae VBR B VE 


which is equation (1), the desired result. 


B. BASS'S METHOD OF DETERMINING A STABLE INITIAL GUESS 

Kleinman's method requires an initial guess of the feedback gain 
that yields a stable closed-loop system control matrix. Kleinman 
remarks [Ref. 2] that, if the system of equation (2) is completely 
controllable, then the desired initial guess will always exist. A 
method that could be programmed for a general, controllable system to 
yield this correct initial guess was sought. 

Bass [Ref. 6] presented, but did not publish, a general method for 
determining a stable initial guess for a completely controllable system 
in 1961. The method was published in a paper by Wonham and Cashman 
[Ref. 7] and a proof can be found in a subsequent paper by Bass [Ref.8]. 


Given the controllabe matrix pair (A,B), the matrix: 
Ay =A - BL, 


] 


will be stable when L, = B'X , where X is the (unique) positive definite 


solution to the linear matrix equation: 
- (A + BI) X - X(A + BI)' + 2BB' = 0 (24) 
where 8 is defined by: 


8 > | {Al 





where ||-|| is the Euclidean norm. According to Wonham [Ref. 7] the 
results are also valid if: 
MAX 
Bae Ne j Ly | a5 5 | 


where as; are the elements of A. 


C. SOLVING THE MATRIX EQUATION Y'X + XY =Z 


Equations (16) and (24) are in the form of the Lyapunov equation: 
Vik tek 2 (25) 


The methods found in the literature for solving this equation fall into 
two general categories: series solutions and simultaneous linear 
equation solutions. 

In the series solutions, X is found from the sum of a converging 
matrix series, 1.e., Ref. 9. For the series to converge the matrix Y 
must be a stable matrix. In solving equation (16) this condition is 
met ; AL is stable by the definition of a bounded cost matrix. In 
general, however, the matrix (A + g1)' of equation (24) is not stable, 
and the series method fails to solve equation (24) properly. 

Since the unique solution, X, is symmetric, equation (25) represents 
n(n + 1)/2 unknowns. The second category of solutions expands equation 
(25) into a set of n(n + 1)/2 simultaneous linear equations. An 
economical way of recursively expanding equation (25) was given by 
Bingulac [Ref. 10], using an integer coefficient matrix to expand the 
equation. This is the method used to expand equations (16) and (24) 
to the form Ax = B, which can be solved by a variety of simultaneous 


equation solvers. 


14 





D. ALTERNATE ANALYTICAL METHODS 

There are several alternate methods found in the literature for 
solving the steady-state matrix Riccati equation (1). 

A method developed by Bass [Ref. 11] obtains the solution from a 
2n-dimensional Hamiltonian, H, and the terms of the polynomial expansion 
of the stable roots of H. This scheme is considered too sensitive to 
finite numerical computations to be of practical use. 

MacFarlane [Ref. 12] shows that the solution can be obtained from 
the eigenvectors corresponding to the unstable eigenvalues of a 
Similar Hamiltonian. The scheme requires that the system have distinct 
eigenvalues and that the corresponding eigenvectors be determined (this 
is difficult for high order systems). 

Blackburn [Ref. 13] programmed a method based on a Newton-Raphson 
iterative solution for simultaneous nonlinear equations. This scheme 
requires an initial guess that yields a stable closed-loop system con- 
trol matrix. The user must determine this initial guess so that it is 
close enough to the final solution for the local Newton-Raphson method 
to converge. 

A fourth method integrates the full, nonsteady-state Riccati 
equation (1) backwards in time, froma set of zero initial conditions, 
until each element of the p matrix reaches a satisfactory, small value. 
For systems of even moderate order, this method is prohibitive with 


respect to computation time. 


Is 





TTT. COMPUTATIONAL METHOD OF SOLUTION 


A. SUBROUTINE RICATS 
1. General 

Subroutine RICATS is the subroutine that the user will call to 
solve equation (1). The program language is FORTRAN IV for the Operat- 
ing System / 360 which is compatible with, and encompasses USASA 
FORTRAN. The calling arguments, in the required sequence, are explained 
in the comment cards at the beginning of the subroutine. It should be 
noted that IA = n and JB =m. Basically the subroutine iterates on 
equation (16), using equation (24) to calculate the initial guess, until 
the solution converges. 

For first order systems (n = 1), the nonlinear equation (1) 
becomes a quadratic equation. For these systems, subroutine RICATS 
solves the resulting quadratic and returns the largest root in the 
first element of the P matrix . 

2. The System Controllability Check 

An analytical requirement of Bass‘'s method of determining the 
initial guess is that the matrix pair (A,B) be completely controllable. 
To check this requirement the augmented matrix G of equation (6) is 
formed and subroutine GMRANK is called to determine the rank of G. If 
the rank equals n, execution continues, otherwise subroutine RICATS 
returns with IER = 2. 

3. The Initial Guess Stability Check 

Kleinman's iteration scheme requires that the eigenvalues of the 

closed-loop system control matrix of the initial guess have negative 


real parts in order for the initial cost matrix UR to be bounded. 
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In subroutine RICATS, when the stability check is requested, the matrix 
A - BL, is copied into a work matrix and then passed to the IBM/SSP 
subroutines HSBG and ATEIG. The eigenvalues computed by ATEIG are then 
individually tested to be algebraically less than minus the absolute 
value of EIGMAX, where EIGMAX is set by the user. If any eigenvalue 
fails the check RICATS returns with IER = 3. 

4. The Solution Positive Definite Check 

Kleinman's iteration scheme is supposed to converge to a positive 

definite solution. The iterations can converge to a nonpositive definite 
solution if all of Kleinman's theorem requirements are not numerically 
met. Therefore, for the user's convenience, a check to verify that the 
solution is positive definite is included. In subroutine RICATS, when 
the positive definite check is requested, the solution is factored by 
the Cholesky square-root method. A nonpositive term on the diagonal 
of the factorized matrix constitutes a nonpositive definite solution 
and RICATS returns with IER = 4 + KL, where KL = 0 is for a converged 
solution and KL = 1 is for NIRY iterations without convergence. 

9. the Steps of Subroutine RICATS 

Calling subroutine RICATS causes the following sequential 

Steps to be executed. 

1. Check input parameters IA, JB, IER, NTRY to insure they are 
within proper bounds. If check fails IER = 6, NN = O and RICATS returns. 
2. If the user requests, check the controllability of the 

inputed system. If check fails IER = 2, NN = O and RICATS returns. 
3. Set NN = 0. 
ieee li this is a first order system (n = 1), go to step 30. 
] 


Seeeiromi ec == BR B’. 
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Seehof = 2Bb. 
7. Form P = (A + gI)' using equation (24) to define g where 


the "greater than" magnitude is provided by the variable FIX or: 


peri | pass |. 

8. Solve (A + BI)X + X(A + BI)' = 2BB'. 

9. If subroutine SIMQ, through subroutine MLIAPS, returns a 
nonzero error flag, NN = 1. 

10. Form L, = Bx”, 

ll. Form P =A - BL, 

12. If the user requests, check the stability of the system 
matrix of the initial guess. If the check fails, IER = 3 and RICATS 
returns. 

13. Form F = - Q - L RL, 

14. Solve (A - BL) ‘Vy a; V5 (A - BL) =-Q- LGRL,- 

15. If subroutine SIMQ, through subroutine MLIAPS, returns a 
nonzero error flag, NN = NN + 1. 

omemoeuoke = 1. KL = 0. 

lf 2 born P= EV, . 
18. Form F =-Q + V_P. 

[oe eOrmer = iA + P. 
+ V 


20. Solve (A + EV.) 'V (A + EV.) So ae eli 


k+] k+] 
21. If subroutine SIMQ, through subroutine MLIAPS, returns a 

nonzero error flag, NN = NN +1. If NN > (n + 1)/2, go to step 29. 
22. Check each element of V,,, by ABS ( (vy, - Vea /Vy ) 

< TOLER. If all elements of Vey Hasseunmise test, go to step 2/. 


C3 FOr VE = Veg 





24. If k > NTIRY; go to step 26. 

Cope oct = ktl: go to step 17. 

26. Set KL = 1. 

27. If the user requests, check to see if Vegy 1S a positive 
definite matrix. If check fails IER = 4 + KL and RICATS returns. 

28. Set JER = KL and RICATS returns. 

29. Set IER = 7 and RICATS returns. 


30. If the user requests the controllability check and B(1) is 


- 


zero, IER = 2, NN = O and RICATS returns. 


ot. set P 


AR/BB + #(AR/BB)* + C'CR/BB 


32. If the user requests the pesitive definite check and P <s 
O.1E-35, IER = 4, NN = O and RICATS returns. 
33. Set IER = 0 and RICATS returns. 


B. SUBROUTINE MLIAPS 

Subroutine MLIAPS is an auxiliary routine used by subroutine RICATS. 
The subroutine expands the equationin steps 8, 14 and 20 of subroutine 
RICATS into n(n + 1)/2 simultaneous linear equations of the form Ax = B. 
The method used was presented by Bingulac [Ref. 10] in 1970. These 
n(n = 1)/2 simultaneous linear equations are then solved by subroutine 
SIMQ. Upon return from subroutine SIMQ, MLIAPS immediately returns to 
Subroutine RICATS and any error codes returned by SIMQ are passed, 


unchanged, to RICATS. 


C. SUBROUTINE GMRANK 
Subroutine GMRANK is a general subroutine for determining the mini- 
mum row or column rank of an arbitary real matrix. The method used is 


Simple row and column interchanges for maximum pivoting, with successive 


he 





reduction on the remaining matrix elements [Ref. 14]. When the absolute 
value of a pivot element is less than 0.1E-35, or when the final pivot 
element has been determined, the subroutine returns the integer K, 

where K is the number of successfully determined nonzero pivot elements 


or, equivalently, the rank of the inputed matrix. 


20 





IV. OPERATIONAL ASPECTS OF USING SUBROUTINE RICATS 


A. CONVERGENCE PROPERTIES 

Kleinman suggests that equation (16) will converge to a satisfactory 
solution within seven to ten iterations, using an initial guess generated 
by hand or from a previous solution. Subroutine RICATS, using a correct 
set of input parameters and Bass's method of generating an initial guess, 
converged within forty iterations for every system tested. 

To test subroutine RICATS, over seventy-five systems were run and a 
Satisfactory solution was returned for each one. As a test of how 
accurate the solutions were, the average absolute value of the matrix p 


of equation (1) was determined. For the satisfactory solutions the 


average values were of the order of 107°. The standard input parameters 
were: 
NTRY = 50 
TOLER = 0Q.1E-03 
FIX = 1.1 
EIGMAX = 0.001 
IER = 0 


and the usual error flag returned was IER = 0 (see the discussion of 


the parameter IER for high order systems). 


B. INPUT PARAMETERS 

Through its input parameters subroutine RICATS was written to be as 
flexible as possible for the user. For any system, the five parameters 
NTRY, TOLER, FIX, EIGMAX, IER influence the final returned solution. 
It is hoped that the user can tailor these input parameters to meet 
his particular needs. 

In the following discussion, “high" or "higher order systems,” 


refer to systems of the eighth order and above (n 2 8). 
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le TUNIS § 
NTRY is the maximum number of iterations the user desires to 
be attempted, before returning to the user's program without a converged 
solution. A recommended value is fifty and the user is usually assured 
that a solution that is going to converge, will converge within fifty 
iterations. It should be noted that some initial guesses generated may 
yield marginally stable system control matrices, and these solutions will 
require a larger number of iterations to converge. From experience, one 
hundred and fifty iterations was considered to be a sufficient practical 
upper limit. 
2. TOLER 
TOLER is the maximum percentage difference, between each element 
of the solution, on successive iterations for the convergence criterion 
to be satisfied. The accuracy of single-precision computations is about 
five significant digits, so decreasing TOLER beyond 0.0001 will not 
result in an appreciably more accurate solution. Decreasing the magni- 
tude of TOLER (up to this limit) will tend to increase the accuracy 
and the number of iterations required for the desired solution. 
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Che FLA, 

FIX is the constant used to insure the magnitude of 8 in equa- 
tion (24) is greater than the Euclidian norm of the A matrix. 
Analytically then, FIX should be greater than one and the suggested 
value of 1.1 worked well for the majority of systems. However, for 
some systems equation (24) generates a singular initial guess, evidenced 
by an underflow error message from subroutine MINV. These singular 


initial guesses can still yield quite satisfactory solutions within 


five or six iterations. However, the user can, by varying FIX through 


Za 





the suggested range, remove this underflow from his execution. From 
practical experience, the range of FIX is from 0.1 up to about three 
Om four. 
4. EIGMAX 

Once the initial guess L is calculated, the user has the option 
through IER of verifying that the associated control matrix is stable. 
If the user asks for the stability check, the eigenvalues of the control 
matrix are determined. Each must have its real part more negative than 
minus the absolute value of EIGMAX. From the discussion on cost matrices, 
theoretically EIGMAX could be set to zero. Numerically though, an 
eigenvalue that 1s very close to zero can induce numerical instability 
in the iteration scheme. From experience, the suggested value is 0.001. 
The user can also use FIX to modify the control matrix of the initial 
guess in an attempt to make the real parts of the eigenvalues negative 
enough to pass the stability check. 

oe «LER 

The most informative parameter of the group is IER. On return 
from RICATS, IER can inform the user of the validity of the solution. 
When speicfying IER as in input parameter, the user can, by bypassing 
various combinations of the three checks available in RICATS, overcome 
some system deficiencies, save execution time, or obviate decks for 
the subroutines GMRANK, ATEIG, HSBG (see note 6 in comment cards at the 
beginning of subroutine RICATS). By bypassing any or all of the three 
checks (if the user is fairly certain that the circumvented checks 
would have been passed ) the user can save execution time, although 


the savings are neglibile except for high order systems. 
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Care must be exercised when not checking controllability and 
or stability. For systems that are uncontrollable, neglecting the 
controllability check may lead to a convergent solution, if the stability 
check were passed. Such a solution could be positive definite and yield 
a stable final-solution control matrix, but the user should keep in 
mind that there may be modes of the system that cannot be controlled. 
For the case of unstable control matrices corresponding to the initial 
guess, the stability check should be bypassed only after FIX has been 
varied through the suggested range of values without success. The 
danger in attempting to generate a solution for a system that would fail 
the stability check is that successive iterations are no longer bounded 
by a lower positive [semi] definite iteration. The iterations will 
probable not converge to a satisfactory answer. This is the most 
dangerous of the three checks to bypass, by far. 

The positive definite check of the solution is included as a 
convenience for the user, to verify that the final solution is indeed 
positive definite. A note of caution should be sounded when solving 
high order systems. The positive definite check is subject to numerical 
problems when evaluating the high order solutions, and thus a solution 
will be flagged as not being positive definite, when for all practical 
purposes it is. Refinements in the solution, from introducing iterative 
routines for solving equation (16) (see subroutine SIMQ below), have 
Shown that the difference between passing and failing the check can 
depend solely on numerically insignificant digits in a few elements of 
the solution. Therefore, to give confidence to a solution that has been 
flagced as not being positive definite, the user can: (1) determine 


the eigenvalues of the final solution, closed-loop system control 
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matrix and check for all negative real parts; or (2) calculate the value 
of the optimal quadratic cost function for an arbitrary set of initial 
conditions using equation (15), and check for a positive, finite result. 
The positive definite check will yield the same result as the IBM/SSP 


Subroutine MFSD. 


Wee HIGH ORDER SYSTEMS 
As the order Of the system increases, the problems due to finite 
numerical calculations increase considerably. One means of trying to 
maintain sufficient accuracy is to have a double-precision version of 
subroutine RICATS. Since this would nearly double the storage require- 
ments (prohibitive for systems of order higher than twenty) it was felt 
that this was not a feasible means of achieving the goal. The suggested 
procedure for systems of higher order is for the user to introduce his 
own dummy subroutines MINV and SIMQ. Care must be exercised when writing 
your own subroutines. The variables passed to and returned from your 
Subroutine must correspond exactly to the variables as handled by the 
Subroutine you are replacing. 
1. Subroutine MINV 

Using a common statement from the user's main program to provide 
the additional storage required, the user can write a routine to convert 
the matrix to be inverted to double-precision storage, invert the matrix 
in double-precision, then convert the inverted matrix back to the 
Single-precision mode. When this is done, the inverted matrix is 
passed back to the subroutine RICATS as if the IBM/SSP subroutine MINV 
had done the inverting. A logical way to invert the matrix in the double- 


precision mode is to use the IBM/SSP subroutine DMINV. It might appear 
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that this type of routine would make no noticeable change in the 
execution of subroutine RICATS. However, for some systems tested, 
particularly those with singular initial guesses (FIX = 1.1), the number 
of iterations required for solution was halved. 
2. Subroutine SIMQ 

Again using a common statement, recognizing that work areas can 
be declared single-precision in one subroutine and double-precision in 
another, the user can write his own simultaneous linear-equation-solving 
routine. The IBM/SSP subroutine SIMQ has been found to yield satis- 
factory results for a maximum of about forty-five equations (n = 9). 
For systems of higher order, the solutions did converge, but they were 
usually accompanied by an error flag indicating that they were not 
positive definite solutions. However, iterative refinement of the 
solution to equation (16) at each iteration can yield error-free 
results. The user can either write his own routine for the iterative 
refinement, or pass the set of equations to the IBM/SSP subroutines 
FACTR and RSLMC. It is suggested that the common statement also contain 
a relative tolerance parameter, not necessarily the same as TOLER 
inputed to RICATS, in either name or magnitude. As was previously 
mentioned, this subroutine technique was used to remove error flags 


from the solutions returned by subroutine RICATS. 
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